O cancelamento de assinatura é uma pedra no sapato de empresas de telecomunicações e streaming, por exemplo, pois traz instabilidade financeira e prejudica o planejamento estratégico da companhia. Por isso, fazer um acompanhamento da taxa de cancelamento ou Churn Rate é fundamental para manter o negócio vivo e de pé. Mas além de monitorar o indicador, é possível prever quando um cliente deve cancelar sua assinatura? Sim, é possível e é isso que vamos mostrar logo mais. Seja bem vindo e te desejo uma boa leitura!
Os dados que vamos utilizar para aplicar aprendizado de máquina estão aqui.
Esta análise tem por finalidade prever o cancelamento de assinatura por parte do cliente, aplicando um modelo de aprendizado de máquina.
library(DescTools)
library(tidyverse)
library(tidymodels)
library(reactable)
library(janitor)
library(highcharter)
library(sparkline)
library(corrplot)
library(patchwork)
library(rmarkdown)
library(knitr)
library(kableExtra)
library(upstartr)
library(tictoc)
library(vip)
theme_set(theme_light(18))
dados = tibble(read_csv("C:/Users/icaro/Documents/Meu_Portfolio/Churn_Telecom/WA_Fn-UseC_-Telco-Customer-Churn.csv"))
O que temos dentro desse conjunto?
paged_table(dados)
Então temos um dataset com 7043 registros/observações e 21 atributos ou melhor, variáveis.
Vamos explorar nosso dataset, começando por conhecer nossas variáveis.
# contagem
vars =
tibble(
tipos = c("quantitativa", "qualitativa"),
colunas =
c(
dados %>% select(where(is.numeric)) %>% ncol(),
dados %>% select(where(is.character)) %>% ncol()
))
# porcentagem
vars =
vars %>%
mutate(
porcent = round(colunas / sum(colunas), digits = 2)
)
# rótulos
rotulo = paste(str_to_title(vars$tipos), " (", vars$porcent * 100, "%", ")", sep = "")
# gráfico
pie(
x = vars$colunas,
labels = rotulo,
col = c("skyblue2", "skyblue4"),
border = "white"
)
title(
main = "Tipo de Variáveis",
cex.main = 2,
font.main = 4,
col.main = "skyblue4"
)
Temos que em sua maior parte, nosso conjunto de dados é formado por variável qualitativa nominal.
Antes de iniciarmos a análise exploratória, vamos realizar alguns ajustes no dataset.
Excluir o campo de ID e padronizar o nome das variáveis.
Transformando o campo Senior Citizen em string.
# padronizando os nomes da vars, excluir a coluna com o ID
dados =
dados %>% select(-customerID) %>% clean_names()
# formatando o campo senior citizen
dados$senior_citizen = ifelse(dados$senior_citizen == 0, 'No', 'Yes')
Agora sim, vamos realizar a construção de um sumário separando entre variáveis numérica e categórica.
# Sumário - Variáveis Categóricas
sumario_character =
dados %>%
select(where(is.character)) %>%
as.list() %>%
enframe(name = "variavel", value = "valores") %>%
mutate(
n_missing = map_dbl(valores, ~sum(is.na(.x))),
complete_rate = 1 - n_missing/map_dbl(valores, ~length(.x)),
min = map_dbl(valores, ~min(str_length(.x))),
max = map_dbl(valores, ~max(str_length(.x))),
# empty = map_dbl(valores, ~sum(str_detect("^$", .x))),
n_unique = map_dbl(valores, ~n_distinct(.x)),
whitespace = map_dbl(valores, ~sum(str_detect("^[:blank:]+$", .x)))
) %>%
arrange(variavel != "churn")
sumario_character %>%
reactable(
wrap = FALSE,
resizable = TRUE,
fullWidth = TRUE,
defaultColDef = colDef(width = 90),
columns = list(
valores = colDef(show = FALSE),
variavel = colDef("Variável", minWidth = 230, width = 230)
),
details = function(index) {
variavel_chr <- sumario_character[index, "variavel", drop = TRUE]
dados %>%
tabyl(!!sym(variavel_chr)) %>% arrange(desc(n)) %>%
reactable(columns = list(percent = colDef("%", format = colFormat(percent = TRUE, digits = 1))), width = 500)
}
)
Agora para as variáveis numéricas.
# Sumário - Variáveis Numéricas
sumario_numeric <- dados %>%
select(where(is.numeric)) %>%
as.list() %>%
enframe(name = "variavel", value = "valores") %>%
mutate(
n_missing = map_dbl(valores, ~sum(is.na(.x))),
complete_rate = 1 - n_missing/map_dbl(valores, ~length(.x)),
mean = map_dbl(valores, ~mean(.x, na.rm = T)),
sd = map_dbl(valores, ~sd(.x, na.rm = T)),
min = map_dbl(valores, ~min(.x, na.rm = T)),
median = map_dbl(valores, ~median(.x, na.rm = T)),
max = map_dbl(valores, ~max(.x, na.rm = T))
) %>%
mutate(across(where(is.numeric), round, digits = 1))
sumario_numeric %>%
reactable(
wrap = FALSE,
resizable = TRUE,
fullWidth = TRUE,
defaultColDef = colDef(width = 90),
columns = list(
valores = colDef(cell = function(values) {sparkline(table(cut(values, min(n_distinct(values), 15))), type = "bar")}),
variavel = colDef("Variável", minWidth = 230, width = 230)
),
details = function(index) {
hchart(sumario_numeric[index, "valores"][[1]][[1]])
}
)
# Correlações entre as numéricas
dados %>%
select(where(is.numeric)) %>%
filter(!is.na(dados$total_charges)) %>%
cor() %>%
corrplot(method = "color",
order = "hclust",
type = "upper",
addCoef.col = "black",
tl.col = "black",
cl.align.text = "l")
# Sankey das explicativas categóricas
dados %>%
select(where(is.character), -churn) %>%
data_to_sankey() %>%
hchart("sankey")
# base transformada
churn_numeric_long =
dados %>%
select(where(is.numeric), churn) %>%
pivot_longer(-churn, names_to = "variavel", values_to = "valor")
# plot base
plot_base =
churn_numeric_long %>%
ggplot(aes(x = valor)) +
facet_wrap(~variavel, scale = "free", ncol = 1) +
theme(axis.text = element_blank(),
axis.title = element_blank())
# Histogramas
p1 <- plot_base + geom_histogram(aes(fill = churn), position = "identity")
# Densidades
p2 <- plot_base + geom_density(aes(fill = churn), alpha = 0.4, show.legend = FALSE)
# Boxplots
p3 <- plot_base + geom_boxplot(aes(fill = churn), show.legend = FALSE)
# KS
p4 <- plot_base + stat_ecdf(aes(colour = churn), show.legend = FALSE)
# plot all
(p4 + p3 + p2 + p1) + plot_layout(nrow = 1)
churn_character_long =
dados %>%
select(where(is.character), churn) %>%
pivot_longer(-churn, names_to = "variavel", values_to = "valor")
# bar plot base
bar_plot =
churn_character_long %>%
count(variavel, valor, churn) %>%
group_by(variavel, valor) %>%
mutate(
p_churn = sum(n[churn == "sim"], na.rm = TRUE)/sum(n, na.rm = TRUE)
) %>%
group_by(variavel) %>%
mutate(
valor = reorder_within(valor, p_churn, variavel)
) %>%
ggplot(aes(y = valor, x = n, fill = churn)) +
facet_wrap(~variavel, scales = "free", ncol = 1) +
scale_y_reordered() +
theme(axis.text.x = element_blank(),
axis.title = element_blank())
# Position stack
p1 = bar_plot + geom_col(show.legend = FALSE, position = "stack")
# Position fill
p2 = bar_plot + geom_col(show.legend = TRUE, position = "fill")
p1 + p2
Vimos que há variáveis categóricas que precisam de alteração para fator, além agrupar os casos de clientes que não possuem serviço de internet ao grupo dos ‘NO’, pois representam a “mesma” situação.
# tranformando em fator as variaveis categoricas
dados_transf =
dados %>%
mutate(
across(
c(gender:dependents,
phone_service, multiple_lines,
online_security:streaming_movies,
paperless_billing, churn),
factor),
multiple_lines =
ifelse(multiple_lines == 'No phone service', 'No',
ifelse(multiple_lines == 'No', 'No', 'Yes')),
online_security =
ifelse(online_security == 'No internet service', 'No',
ifelse(online_security == 'No', 'No', 'Yes')),
online_backup =
ifelse(online_backup == 'No internet service', 'No',
ifelse(online_backup == 'No', 'No', 'Yes')),
device_protection =
ifelse(device_protection == 'No internet service', 'No',
ifelse(device_protection == 'No', 'No', 'Yes')),
tech_support =
ifelse(tech_support == 'No internet service', 'No',
ifelse(tech_support == 'No', 'No', 'Yes')),
streaming_tv =
ifelse(streaming_tv == 'No internet service', 'No',
ifelse(streaming_tv == 'No', 'No', 'Yes')),
streaming_movies =
ifelse(streaming_movies == 'No internet service', 'No',
ifelse(streaming_movies == 'No', 'No', 'Yes'))
)
Vamos iniciar o processo de modelagem dos dados com a separação da base em teste e treino.
# marcando a base em teste e treino
set.seed(1232)
churn_split =
initial_split(data = dados_transf, strata = churn, prop = 0.75)
# separando a base em treino
training_df = training(churn_split)
# separando a base em teste
testing_df = testing(churn_split)
# avaliação dos modelos com cross-validation
folds =
vfold_cv(training_df, v = 10, strata = churn)
Na receita informamos as etapas que precisamos realizar para preparar a base de treino do modelo.
No nosso caso, a receita inclui:
Correção do balanceamento de classes.
Exclusão de variáveis com alto coeficiente de correlação.
Normalização das variáveis numéricas.
Substituição de valores faltantes pela mediana em campos númericos.
Transformação das variáveis qualitativas nominais em dummy.
# preparando os dados
churn_rec =
recipe(churn ~ ., data = training_df) %>%
themis::step_downsample(churn) %>%
step_corr(all_numeric_predictors()) %>% # correlação entre explicativas
step_normalize(all_numeric()) %>% # normalizar vars
step_rm(total_charges) %>% # Tratamento de missings
step_impute_median(all_numeric()) %>% # Substituir os missings com a mediana
step_dummy(all_nominal(), -churn) # Transformando as variaveis nominais em dummy
# base processada
baked = bake(prep(churn_rec), new_data = NULL)
paged_table(baked)
# modelo de randow forest (árvore aleatória)
rf_mod =
rand_forest(
trees = 300,
mtry = tune(),
min_n = tune()
) %>%
set_mode('classification') %>%
set_engine('ranger', importance = 'impurity')
rf_workflow =
workflow() %>%
add_recipe(churn_rec) %>%
add_model(rf_mod)
set.seed(1)
tic("modelo rf")
tunagem =
tune_grid(
rf_workflow,
resamples = folds,
grid = 10,
metrics = metric_set(roc_auc, precision, accuracy, f_meas),
control = control_grid(verbose = TRUE, allow_par = FALSE)
)
toc()
## modelo rf: 79.81 sec elapsed
# tabela
show_best(tunagem, "roc_auc") %>%
kable(format = 'html', align = 'l', digits = 3) %>%
kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover", "condensed", "responsive"))
| mtry | min_n | .metric | .estimator | mean | n | std_err | .config |
|---|---|---|---|---|---|---|---|
| 8 | 31 | roc_auc | binary | 0.848 | 10 | 0.005 | Preprocessor1_Model02 |
| 13 | 36 | roc_auc | binary | 0.848 | 10 | 0.005 | Preprocessor1_Model01 |
| 22 | 38 | roc_auc | binary | 0.845 | 10 | 0.005 | Preprocessor1_Model10 |
| 4 | 8 | roc_auc | binary | 0.844 | 10 | 0.006 | Preprocessor1_Model09 |
| 17 | 26 | roc_auc | binary | 0.844 | 10 | 0.005 | Preprocessor1_Model04 |
Atualizando o workflow com os hiperparamêtros encontrados.
# atualizando wf
wf =
rf_workflow %>%
finalize_workflow(select_best(tunagem, "roc_auc"))
# last fit
ajuste_final =
last_fit(
wf,
churn_split,
metrics =
metric_set(accuracy, roc_auc, f_meas, specificity, precision, recall)
)
# métricas de desempenho
collect_metrics(ajuste_final) %>%
kable(format = 'html', digits = 3, align = 'l') %>%
kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover", "condensed", "responsive"))
| .metric | .estimator | .estimate | .config |
|---|---|---|---|
| accuracy | binary | 0.740 | Preprocessor1_Model1 |
| f_meas | binary | 0.805 | Preprocessor1_Model1 |
| specificity | binary | 0.769 | Preprocessor1_Model1 |
| precision | binary | 0.897 | Preprocessor1_Model1 |
| recall | binary | 0.730 | Preprocessor1_Model1 |
| roc_auc | binary | 0.825 | Preprocessor1_Model1 |
# predicoes
pred_base_teste = collect_predictions(ajuste_final)
# curva roc
roc =
pred_base_teste %>%
roc_curve(churn, .pred_No) %>%
autoplot()
# curva de lift
lift =
pred_base_teste %>%
lift_curve(churn, .pred_No) %>%
autoplot()
# KS
ks =
pred_base_teste %>%
ggplot(aes(x = .pred_Yes, colour = churn)) +
stat_ecdf(show.legend = FALSE)
# distribuicao
dist =
pred_base_teste %>%
ggplot(aes(x = .pred_Yes, fill = churn)) +
geom_density() +
theme(axis.title = element_blank())
(roc + lift)/(ks + dist)
A conclusão que tiramos daqui é que o modelo apresentou razoável acurácia, separou bem os grupos, mas há espaço para melhorar, pois existe um ponto central no gráfico de densidade que mostra que o modelo não conseguiu acertar.
Ajuste final com a base inteira.
modelo_final = fit(wf, dados_transf)
Avaliação das contribuição das variáveis para o modelo.
vip(modelo_final$fit$fit)
Podemos salvar o modelo para realizar as predições futuramente, sem a necessidade de rodar o script sempre que precisar.
saveRDS(modelo_final, file = "modelo_final.rds")
Vamos confrontar as previsões do modelo com a realidade.
predizer =
tibble(
cbind(
dados_transf$churn,
predict(modelo_final, new_data = dados_transf, type = 'prob')
)
)
predizer =
predizer %>%
rename(churn_observado = `dados_transf$churn`) %>%
mutate(churn_esperado = ifelse(.pred_No > .pred_Yes, 'No', 'Yes'),
.pred_No = round(.pred_No, digits = 3),
.pred_Yes = round(.pred_Yes, digits = 3))
# matriz
predizer %>%
tabyl(churn_observado, churn_esperado) %>%
adorn_percentages(denominator = "all") %>%
adorn_pct_formatting() %>%
adorn_title(row_name = "Observado", col_name = "Esperado") %>%
kable(format = "html", align = 'c', caption = "Matriz de confusão") %>%
kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover", "condensed", "responsive"))
| Esperado | ||
|---|---|---|
| Observado | No | Yes |
| No | 56.9% | 16.6% |
| Yes | 3.5% | 23.0% |
Vamos testar a predição do modelo, inserindo os dados de um ‘novo cliente’.
novo_cliente =
data.frame(
gender = 'Female',
senior_citizen = 'Yes',
partner = 'Yes',
dependents = 'Yes',
tenure = 25,
phone_service = 'No',
multiple_lines = 'No',
internet_service = 'DSL',
online_security = 'No',
online_backup = 'No',
device_protection = 'No',
tech_support = 'No',
streaming_tv = 'Yes',
streaming_movies = 'Yes',
contract = 'Month-to-month',
paperless_billing = 'Yes',
payment_method = 'Electronic check',
monthly_charges = 55.55,
total_charges = 290.50
)
predict(modelo_final, new_data = novo_cliente, type = 'prob') %>%
kable(
format = "html",
digits = 3,
align = "c",
caption = "Novo cliente",
col.names = c("Não", "Sim")
) %>%
kable_styling(
full_width = FALSE,
bootstrap_options = c("striped", "hover", "condensed", "responsive")
)
| Não | Sim |
|---|---|
| 0.509 | 0.491 |
Neste exemplo, com a probabilidade maior para não, então assumimos que o cliente não deve cancelar sua assinatura.
Recaptulando o fizemos aqui até.
Análise descritiva da base, explorando as relações entre as variáveis.
Transformação necessária para aplicação do modelo.
Modelagem, com a especificação da f(x) - receita -
workflow até o armazenamento do modelo e avaliação.
Ao final deste estudo, podemos concluir que o modelo conseguiu realizar um trabalho razoável, com precisão de 0,89. No entanto, a matriz de confusão nos mostrou que o modelo errou bastante, onde previu que o cliente cancelaria a assinatura quando na verdade não cancelou.